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Spatial Effects of Delay-induced Stochastic Oscillations in a Multi-Cellular System 


Dmitry Bratsun 1 and Andrey Zakharov 2 

1 Department of Applied Physics, Perm National Research Polytechnical University, 614990, Perm, Russia 
2 Department of Chemical Engineering, Technion - Israel Institute of Technology, 32000, Haifa, Israel 

We explore the joint effect of the intrinsic noise and time delay on the spatial pattern formation 
within a multi-scale mobile lattice model of the epithelium. The protein fluctuations are driven by 
transcription/translation processes in epithelial cells exchanging chemical and mechanical signals 
and are described by a single-gene auto-repressor model with constant delay. Both deterministic 
and stochastic descriptions are given. We found that time delay, noise and spatial signaling can 
result in the protein pattern formation even when deterministic description exhibits no patterns. 
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The small number of reactant molecules involved in 
gene regulation can lead to significant fluctuations in pro¬ 
tein concentrations, and there have been numerous stud¬ 
ies devoted to the influence of such noise at the regulatory 
level since pioneering works in the early 2000s [lH3] . For 
a good review of recent developments, see 00. 

The transcription-translation processes are compound 
multistage reactions involving the sequential assembly of 
long molecules. It can provoke a time lag in gene regu¬ 
lation processes. Until delays are small compared with 
other significant time scales characterizing the genetic 
system, one can safely ignore them in simulations. How¬ 
ever, if the lags become longer than other processes, the 
system has to be considered as non-Markovian, and one 
should account for it in both deterministic and stochas¬ 
tic descriptions. The joint effect of the intrinsic noise and 
time delay on the temporal behavior during gene regula¬ 
tion have been studied first in EH3. We have suggested a 
generalization of the Gillespie algorithm widely used 

to simulate statistically correct trajectories of the state of 
a chemical reaction network that accounts for delay [6J. 
Based on this technique, we showed that quasi-regular 
fluctuations can arise in the stochastic system with delay 
even when its deterministic counterpart exhibits no os¬ 
cillations [71. Since that papers, there have been a lot of 
works developing this line of research (see [iCEHE] for 
recent reviews). They mainly focus on further improving 
the algorithm and on studying the temporal dynamics of 
the different gene systems with delays. 

Several years ago, Lemerle et al. H3J have noticed that 
“space is the final frontier in stochastic simulations of bio¬ 
logical systems”. The problem is that the massive amount 
of spatio-temporal experimental data has been accumu¬ 
lating, but stochastic models of biochemical processes are 
focusing mostly on the temporal dynamics. If in the past 
years spatial stochastic simulations of Markovian pro¬ 
cesses have made considerable progress El HMS1, the 
examples of studies of non-Markovian stochastic systems 
are very rare. For instance, Marquez-Lago et al. PH dis¬ 
cussed how the spatial displacement of molecules can be 
incorporated into purely temporal models through dis¬ 
tributed delays. Danino et al. [18] have given the re¬ 


markable experimental data with the pattern formation 
of delay-induced rhythms in a population of E. Coli , but 
did not provide the stochastic modeling. The theoretical 
difficulties seems to be clear: the generalization of the 
Gillespie algorithm to the case of the spatial dynamics of 
time-delayed processes is still waiting for its author. 

In this Letter, we explore the spatial effects produced 
by the combined action of the intrinsic noise and time 
delay within a multi-cellular system. It is based on the 
multi-scale chemo-mechanical model of the epithelial tis¬ 
sue which first was developed in m and then was applied 
to simulate the carcinoma growth [20] . Using the lattice 
approach allows to avoid the problem of the lack of a re¬ 
liable Gillespie-like algorithm for spatial non-Markovian 
processes and to apply the algorithm from |6]. We use a 
single-gene auto-repressor model with constant delay at 
a single-cell level complemented by the dynamics of the 
transport species penetrating cell membranes (Fig. [l]). 
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FIG. 1: A single-gene autorepressor oscillator with a time- 
delayed negative feedback. Its protein (red circles) acts as a 
positive regulator of the transport protein (blue circles) by 
activating its transcription. The signaling species is diffused 
from one cell to the other within the lattice governed by the 
chemo-mechanical epithelium model described in [ J9I 20L 
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Epithelial tissue is a two-dimensional layer of cells 
covering the surface of an organ or body. The chemo- 
mechanical model includes the calculation of the dynam¬ 
ics of separate cells adapting to environmental stresses. 
The initial hexagonal lattice representing cells becomes 
distorted due to proliferation/intercalation and eventu¬ 
ally incorporates also polygons with a different number 
of vertices. The epithelium is evolved by moving the cell 
nodes. The mechanical force acting on any ith node is de¬ 
fined via the elastic potential energy of the tissue [19j 20] : 

F, = ii^( Ki2+,)(5 " So ) 2 ) +Ff ’ (1) 

cells 




FIG. 2: (a) Neutral curve of the Hopf bifurcation obtained 
within a single-cell deterministic description for e = 0.1, S = 
0.2. (b) Power spectrum obtained in a single-cell stochastic 
simulation below the Hopf bifurcation for A = 20, B = 5, 
t — 6, fc+i — 100, k+a — 200, k- 1 — 1000, k-d — 1000. 


where is the radious vector of ith. node, L and S stand 
for the perimeter and area of a cell respectively. Ff* is an 
uncorrelated zero mean stochastic input. The coefficient 
k defines the action of contractile forces within the cell 
cortex, rj reflects the cell resistance to any changes with 
respect to the reference area So- 

Since the motion is strongly overdamped, the govern¬ 
ing equation for the displacement should be written as 

-^ = F i H(\F i \-F 0 ), (2) 

where H is a Heaviside function, Fq is the threshold force 
below which the node remains immobile. Altogether, 
Eqs. define the mechanics of the tissue. 

We consider a single-gene protein synthesis with neg¬ 
ative auto-regulation (Fig. [I]). This is a popular motif 
in genetic regulatory circuits, and its temporal dynamics 
has been analyzed within both deterministic and stochas¬ 
tic framework [l]. This model is relatively simple yet but 
still maintains a high degree of biological relevance. Its 
generalized version accounting for the effect of time delay 
r has been suggested in mm- Suppose that protein can 
exist both in the form of monomers X and dimers X^. 
The transitions between them with the rates k±d are 

X + X-^>X D , X d -Hx + X. (3) 

We assume that the protein may be degraded with the 
rate B and be produced with the rate A respectively: 

X A 0, D‘ A D‘ +T + X t+T . (4) 

The synthesis occurs at time t + r if the chemical state of 
the promoter site of the X gene at time t is unoccupied 
(Dq). Otherwise (Di), the production is blocked. The 
transitions between the states occur with rates k±\ dur¬ 
ing binding and unbinding of some dimer respectively 

Do + X^D, D, —4 D 0 + X D . (5) 

In order to describe the intercellular signaling, we in¬ 
troduce the transport species T positively regulated by 


the X protein (Fig. [TJ. If the operator site of the T gene 
is occupied (D^) then the T protein may be produced 
immediately with a certain probability At' 


g A T DT(t) 


» T. 


(6) 


If the site is unoccupied, the production of the T protein 
is blocked. Thus, the X monomers act as a positive reg¬ 
ulator of T by activating its transcription, because the 
transitions between the states with the rates k ±2 are 


dJ + x^dJ, 


N —4 Dq 


X. 


( 7 ) 


We assume also that once a signal has come in a certain 
cell, it is converted into the X monomers with the rate 
Bt and the copy number N\ 


X. 


(8) 


Altogether, Eqs. (|3][8]) define the kinetics of gene regula¬ 
tion both at a single cell level and a whole tissue. 

Deterministic description. The main approximation is 
that the reactions of dimerization ^ and binding/ un¬ 
binding ( 5j7 ) are fast in comparison with production/ 
degradation of proteins ®8j). Thus, their dynamics has 
to enter quickly into a local equilibrium and we arrive to 


N dxi A 

^ + £Xi) ~dt ~ 1 + e5x 2 i {t-T) 


+ B T 0i — Bxi, (9) 


dOi NArcrXi 
dt 1 + axi 


-B T 0i+ aLijiOj-Oi), (10) 

jEadj(i) 


where the subscripts refer to cells, x and 0 stand for the 
concentrations of the X and T proteins respectively, a 
is the transfer coefficient, S = k+d/k-d, £ = fc+iA-i, 
a = /c_|_ 2 / k —2 and adj{i ) stands for “adjacent to i-cell”. 
It is assumed that the T protein is transported diffusively 
from one cell to the other, whereas its flux does not de¬ 
pend on the distance between the two cells i and j but 
is proportional to the boundary length Lij. This implies 
that the transport is limited by the transfer though cell 
membranes. The link between sub-cellular and macro¬ 


scopic scales is established through the Eq. (10), since 
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the 9 field is global for the whole tissue (for more details, 
see Supplementary materials). 

The neutral curve for the Hopf bifurcation derived 
within the deterministic approach ( 9fl0 ) at a single-cell 
level is plotted in Fig. & in the plane of dimensionless 
parameters r A and tB. The numerical study reveals the 
common dynamics below and above the bifurcation. 

In order to study the spatial effects, the set of delay 
differential equations has been solved using the 

explicit Euler method, whose stability was warranted by 
a sufficiently small time step. This procedure was syn¬ 
chronized with the simulation of the mechanical evolu¬ 
tion governed by Eqs. ([T][2|. The initial configuration of 
the system is a hexagonal lattice comprising 1560 cells 
with random phase distribution. The tissue as a whole 
has the form of a stripe. The typical values of the pa¬ 
rameters governing the tissue mechanics hereinafter are 
as follows: k = 1.0, r] = 1.0, So = 3>/3/2, F 0 = 0.02. 
Fig. [3^i presents the results of numerical simulation with 
parameters taken above the Hopf bifurcation curve (the 
upper black square in Fig. [2^i) . The X protein patterns 
are shown for two consecutive moments of time. The non¬ 
linear dynamics includes the slow development of spiral 
traveling wave pattern which arises against the synchro¬ 
nized field oscillating in the background. The oscillation 
period is approximately equal to the triple delay time. 

Stochastic description. In order to describe the spatial 
stochastic effects, we use a hybrid model, which is con¬ 
structed as follows. The dynamics of the proteins in each 
cell has been obtained by performing direct stochastic 
simulations of the reactions ([3][8| using the modified ver¬ 
sion of the Gillespie algorithm j6] . The signaling between 
cells is still organized as diffusive transport from one cell 
to the other according to finite-difference formula: 


/S.t 


= Tj + 


At a Lij(Tj — Tf) 

jeadj(i) 


(ii) 


where [...] stands for the integer part of the expression. 
The time step in (11) is equal to the time step for the in¬ 
tegration of (plf 


At = 0.05. Since the typical time step 
stochastic system is much less (A t st = 0.00001 — 0.001), 
one needs to dock the numerical schemes for the mechan¬ 
ical evolution of the tissue and stochastic fluctuations of 
the protein (see Supplementary materials). 

The power spectrum of time-delayed stochastic signal 
generated on a single-cell level far below the Hopf bi¬ 
furcation is shown in Fig. It. It has a remarkable row 
of peaks indicating that some periodicity in the signal 
should exist. The first peak corresponds to the frequency 
cj* = 7r/r ~ 0.524 (r = 6). One can observe also the 
strong harmonics of the fundamental frequency given by 
uj n = (2 n — 1)cj*. Thus, the time delay coupled with the 
noise can cause a system to be oscillatory even when its 
deterministic description predicts no oscillations. 



FIG. 3: Evolution of the X protein pattern in the epithelium 
formed of more than 1500 cells within (a) deterministic and 
(b) stochastic description, with parameters A = 500, At = 
500, B = 5, r = 6, N = 1, a = 0.05, k d = 200, k- t = 1000, 
An ,2 = 100 taken above deterministic Hopf bifurcation curve 
(see Fig. |2^t). The frames correspond to times t = 260, 340. 


Let us consider now two examples of spatial stochas¬ 
tic simulations. Fig. presents the stochastic pattern 
formed by the X monomers in the tissue with the same 
parameter values as in Fig. UK We found that nonlinear 
dynamics of spatially extended system consists of two dis¬ 
tinct oscillatory modes, just like it was in the determinis¬ 
tic case. One is a quasi-standing wave pattern oscillating 
with u;*. The second oscillatory mode consists of travel¬ 
ing waves which arise from selected initial disturbances. 
In fact, the stochastic pattern looks very similar to its 
deterministic counterpart obtained for the same param¬ 
eters (compare with the upper row in the same figure). 
The wavelength of the structure is found to depend on 
the copy number N whose growth enhances fluctuations 
and diffusive fluxes between cells. 

Consider now the parameters taken far below the Hopf 
bifurcation curve (Fig. [2^). We found that starting with 
random initial conditions, the system fairly quickly falls 
into a fully synchronized mode of oscillations with a com¬ 
mon frequency uj ~ cj* (Fig.[4j[5]). We found also that the 
increase of the copy number N can result in the clustering 
when the cells form two approximately equal communi¬ 
ties, which oscillate in anti-phase. For instance, the nu¬ 
merical simulation with copy number N = 4 has showed 
that the clustering is not observed. In contrast to that, 
this effect manifests itself clearly at N = 8 after a suffi¬ 
ciently long integration (Fig. [4]). In fact, the clustering in 
the system with a large amount of elements exchanging 
chemical signals has become at the center of attention 
of many scientists recently (see, for example, EB). it is 
believed that the clustering is likely to be the reason of 
further cells differentiation in organs. 
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FIG. 4: Stochastic pattern formation of the X protein in the epithelium far below the Hopf bifurcation (see Fig. [2^). The 
frames correspond to times within one oscillatory period 2r. The parameters are the same as in Fig. [3] except A = 20, N = 8. 


Closing remarks. It is known that the noise during 
gene expression conies about in two ways. The inher¬ 
ent stochasticity of transcription/translation generates 
intrinsic noise. The extrinsic noise refers to variation 
in identically-regulated quantities between different cells. 
In this paper, we have focused on spatial effects of in¬ 
trinsic noise. A reliable Gillespie-like algorithm for the 
spatial non-Markovian systems still has not been devel¬ 
oped, but perhaps, this algorithm is not very necessary. 
Since any living matter consists of cells, it is more prac¬ 
tical to use numerical methods based on a computational 
mesh which represents a cellular compartment, such as a 
membrane or the interior of some part of a cell. In this 
paper, we have applied such lattice approach suggesting 
the model which includes both mechanical interactions 
and chemical signal exchanges between cells. 

An important new result of this study is giving insight 
into how the excitation of quasi-regular delay-induced 
fluctuations found in [6], manifests itself in space. We 
show that above the Hopf bifurcation it is observed the 
traveling wave pattern which is similar to that in the de¬ 
terministic case. But more interesting result is found be¬ 
low the Hopf bifurcation where the deterministic system 
is stable: there may be observed both a spatial synchro¬ 
nization of oscillations and clustering of cell community. 

We wish to thank L.M. Pismen for stimulating discus¬ 
sions. The work was supported by the Perm Ministry of 
Education and RFBR (grant 14-01-96022r_ural_a). 
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